#### How Immigration Policies Sustain Authoritarian Regimes: A Case of Saudi Arabia
#### Masaki Matsuo and Shingo Hamanaka (February 2025, revised April 2025)

# Some packages such as "dplyr" and "tidyverse" 
# may cause conflicts with package "list" 
# which is applied in this paper.

rm(list = ls())
load("saudiictreg_2025April.RData")

#### Balance Check (described in the appendices)####
#age
t.test(age~exp1, data=saudiictreg)
#income
t.test(income ~ exp1, data = saudiictreg)
#education
t.test(education ~ exp1, data = saudiictreg)
#Ethnocracy
t.test(Ethnocracy ~ exp1, data = saudiictreg)
#kinship
t.test(Kinship ~ exp1, data = saudiictreg)
#left_behind
t.test(Left_Behind ~ exp1, data = saudiictreg)
#gender
saudiictreg$gender<-ifelse(saudiictreg$gender=="Male",0,1)
t.test(gender ~ exp1, data = saudiictreg)
#Riyad
t.test(Riyad ~ exp1, data = saudiictreg)
#Jeddah
t.test(Jeddah ~ exp1, data = saudiictreg)
#Dammam
t.test(Dammam ~ exp1, data = saudiictreg)

#### Step1: Estimating the support for MbS ####
table(response=saudiictreg$outcome, 
      group=saudiictreg$exp1)
result<-t.test(outcome~exp1, data = saudiictreg)
result$estimate[2]-result$estimate[1]
print(result)
result_ols<-lm(outcome ~ exp1, data = saudiictreg)
summary(result_ols)

#### Step2: Regression Analysis ####
library(list)
full.model<- ictreg(formula=outcome~Ethnocracy 
                    + Kinship 
                    + education 
                    + income,
                    data=saudiictreg,
                    treat="exp1",J=3, 
                    constrained=FALSE, 
                    method="ml")

# Estimating differences between groups
df.ethnocracy <- df.kinship <- saudiictreg
df.ethnocracy[,"Ethnocracy"]<-1;df.ethnocracy[,"Kinship"]<-0
df.kinship[,"Kinship"]<-1;df.kinship[,"Ethnocracy"]<-0
df.left_behind <-saudiictreg
df.left_behind[,c("Ethnocracy","Kinship")]<-0
df.left_behind[,"Left_Behind"]<-1

# Difference between Ethnocracy and Kinship
avg.pred.diff.eth_knsp <- predict(full.model, 
                                  newdata=df.ethnocracy, 
                                  newdata.diff=df.kinship, 
                                  se.fit=TRUE, avg=TRUE, 
                                  interval="confidence")

avg.pred.diff.eth_knsp$fit

# Difference between Ethnocracy and Left_Behind
avg.pred.diff.eth_LB <- predict(full.model, newdata=df.ethnocracy,
                                newdata.diff=df.left_behind,
                                se.fit=TRUE, avg=TRUE,
                                interval="confidence")

avg.pred.diff.eth_LB$fit

# Difference between Kinship and Left_Behind
avg.pred.diff.kinship_LB <- predict(full.model,
                                    newdata=df.kinship,
                                    newdata.diff=df.left_behind,
                                    se.fit=TRUE, avg=TRUE,
                                    interval="confidence")

avg.pred.diff.kinship_LB$fit

# Estimating differences between gender (not in the main text)
full.model.gender <- ictreg(formula=outcome~ gender
                            + education 
                            + income,
                            data=saudiictreg,
                            treat="exp1",J=3, 
                            constrained=FALSE, 
                            method="ml")

df.female <- df.male <- saudiictreg
df.male[,"gender"]<-0
df.female[,"gender"]<-1

# Difference between male and female (not in the main text)
avg.pred.diff.gender <- predict(full.model.gender, 
                                newdata=df.male, 
                                newdata.diff=df.female, 
                                se.fit=TRUE, avg=TRUE, 
                                interval="confidence")

avg.pred.diff.gender$fit

# Estimating differences between generations (not in the main text)
# building variable of young(18 to 24 years old)
saudiictreg$young<-ifelse(saudiictreg$age<=1,1,0)
table(saudiictreg$young)

full.model.young<- ictreg(formula=outcome~ young
                          + education 
                          + income,
                          data=saudiictreg,
                          treat="exp1",J=3, 
                          constrained=FALSE, 
                          method="ml")


# Estimating differences between generations
df.others <- df.young <- saudiictreg
df.young[,"young"]<-1
df.others[,"young"]<-0

# Difference between generations
avg.pred.diff.young <- predict(full.model.young, 
                               newdata=df.young, 
                               newdata.diff=df.others, 
                               se.fit=TRUE, avg=TRUE, 
                               interval="confidence")

avg.pred.diff.young$fit

# Estimating differences between regions (not in the main text)
full.model.region<- ictreg(formula=outcome~Riyad
                           + Jeddah
                           + education 
                           + income,
                           data=saudiictreg,
                           treat="exp1",J=3, 
                           constrained=FALSE, 
                           method="ml")

# Estimating differences between regions
df.Dammam <- df.Jeddah <- df.Riyad <- saudiictreg
df.Riyad[,"Riyad"]<-1;df.Riyad[,"Jeddah"]<-0;df.Riyad[,"Dammam"]<-0
df.Jeddah[,"Jeddah"]<-1;df.Jeddah[,"Riyad"]<-0;df.Jeddah[,"Dammam"]<-0
df.Dammam[,"Dammam"]<-1;df.Dammam[,"Riyad"]<-0;df.Dammam[,"Jeddah"]<-0

# Difference between Riyad and Jeddah
avg.pred.diff.R_J <- predict(full.model.region, 
                             newdata=df.Riyad, 
                             newdata.diff=df.Jeddah, 
                             se.fit=TRUE, avg=TRUE, 
                             interval="confidence")

avg.pred.diff.R_J$fit

# Difference between Riyad and Dammam
avg.pred.diff.R_D <- predict(full.model.region, newdata=df.Riyad,
                             newdata.diff=df.Dammam,
                             se.fit=TRUE, avg=TRUE,
                             interval="confidence")

avg.pred.diff.R_D$fit

# Difference between Jeddah and Dammam
avg.pred.diff.J_D <- predict(full.model.region,
                             newdata=df.Jeddah,
                             newdata.diff=df.Dammam,
                             se.fit=TRUE, avg=TRUE,
                             interval="confidence")

avg.pred.diff.J_D$fit


#### Step3  Estimating preference of three groups for policy ####
library(cjoint)

load("HIPSAR_cjointsaudi_2025April.Rdata")

attr_list <- list()
attr_list[["number"]] <- c("Yes","No")
attr_list[["nationality"]] <- c("Yes","No")
attr_list[["duration"]] <- c("Yes","No")
attr_list[["occupation"]] <- c("Yes","No")
attr_list[["education"]] <- c("Yes","No") 
attr_list[["dependents"]] <- c("Yes","No")
attr_list[["age"]] <- c("Yes","No")
attr_list[["residence"]] <- c("Yes","No")
attr_list[["religion"]] <- c("Yes","No")

# base model
cjoint_design <- makeDesign(type = "constraints",
                            attribute.levels = attr_list)
cjoint_pool_base <- amce(choice ~ number 
                         + nationality 
                         + duration 
                         + occupation 
                         + education 
                         + dependents 
                         + age 
                         + residence 
                         + religion, 
                         data=HIPSAR_cjointsaudi,
                         respondent.id= "ID",
                         cluster = TRUE,
                         design = cjoint_design,
                         na.ignore=TRUE)
summary(cjoint_pool_base)

# Intersection: Ethnocracy
cjoint_ethnocracy <- amce(choice ~ Ethnocracy*(number
                                               + nationality 
                                               + duration 
                                               + occupation 
                                               + education 
                                               + dependents 
                                               + age 
                                               + residence 
                                               + religion),
                          data=HIPSAR_cjointsaudi,
                          respondent.id= "ID",
                          respondent.varying = "Ethnocracy",
                          cluster = TRUE,
                          design = cjoint_design,
                          na.ignore=TRUE)
summary(cjoint_ethnocracy)

# Intersection:Kinship
cjoint_kinship <- amce(choice ~ Kinship*(number
                                         + nationality 
                                         + duration 
                                         + occupation 
                                         + education 
                                         + dependents 
                                         + age 
                                         + residence 
                                         + religion),
                       data=HIPSAR_cjointsaudi,
                       respondent.id= "ID",
                       respondent.varying = "Kinship",
                       cluster = TRUE,
                       design = cjoint_design,
                       na.ignore=TRUE)
summary(cjoint_kinship)

# Intersection: Left-Behind
cjoint_LB <- amce(choice ~ Left_Behind*(number 
                                        + nationality 
                                        + duration 
                                        + occupation 
                                        + education 
                                        + dependents 
                                        + age 
                                        + residence 
                                        + religion),
                  data=HIPSAR_cjointsaudi,
                  respondent.id= "ID",
                  respondent.varying = "Left_Behind",
                  cluster = TRUE,
                  design = cjoint_design,
                  na.ignore=TRUE)
summary(cjoint_LB)

#### drawing figures ####
#### figure 1 ####
#extracting data from the result (avg.pred.diff.eth_knsp) to pass to the modelplot function
library(modelsummary)
library(dplyr)
plotdata.fig1<-avg.pred.diff.eth_knsp$fit
plotdata.fig1<-mutate(plotdata.fig1, 
                      term = c("Ethnocracy","Kinship","Difference"))
plotdata.fig1<-mutate(plotdata.fig1,
                      std.error = avg.pred.diff.eth_knsp$se.fit)
colnames(plotdata.fig1)[1:3]<-c("estimate","conf.low","conf.high")

plotdata.fig1<-list(
  tidy = plotdata.fig1,
  glance = NA
)
class(plotdata.fig1)<-"modelsummary_list"
plotdata.fig1$tidy$model
plotdata.fig1$tidy$model<-c("Ethnocracy", "Kinship", "Difference")

#building data for figure 1

fig1<-modelplot(list("model1" = plotdata.fig1),
                size=0.4)+
  labs(x = NULL)+
  scale_x_continuous(
    limits = c(-0.125, 1.25),
    breaks = seq(-0.5, 1.25, by = 0.25),
    labels = function(x) ifelse(x %% 0.25 == 0, x, "")) +
  geom_vline(xintercept = 0, linetype = 2)+
  theme_bw(base_size = 11)+
  coord_flip()
fig1

#saving fig1 as a TIFF file.
tiff("figure1.tiff", width = 4, height = 4, units = "in", res = 600)
print(fig1)
dev.off()

#### figure 2 ####
#extracting data from the result(avg.pred.diff.eth_LB) to pass to the modelplot function
plotdata.fig2<-avg.pred.diff.eth_LB$fit
plotdata.fig2<-mutate(plotdata.fig2, 
                      term = c("Ethnocracy","Left-Behind","Difference"))
plotdata.fig2<-mutate(plotdata.fig2,
                      std.error = avg.pred.diff.eth_LB$se.fit)
colnames(plotdata.fig2)[1:3]<-c("estimate","conf.low","conf.high")

plotdata.fig2<-list(
  tidy = plotdata.fig2,
  glance = NA
)
class(plotdata.fig2)<-"modelsummary_list"

#building data for figure 2
fig2<-modelplot(list("model2" = plotdata.fig2),
                size=0.4)+
  labs(x = NULL)+
  scale_x_continuous(
    limits = c(-0.5, 1.25),
    breaks = seq(-0.5, 1.25, by = 0.25),
    labels = function(x) ifelse(x %% 0.25 == 0, x, "")) +
  geom_vline(xintercept = 0, linetype = 2)+
  coord_flip()
fig2

#saving fig2 as a TIFF file.
tiff("figure2.tiff", width = 4, height = 4, units = "in", res = 600)
print(fig2)
dev.off()

#### figure 3 ####
#extracting data from the result (avg.pred.diff.kinship_LB) to pass to the modelplot function
plotdata.fig3<-avg.pred.diff.kinship_LB$fit
plotdata.fig3<-mutate(plotdata.fig3, 
                      term = c("Kinship","Left-Behind","Difference"))
plotdata.fig3<-mutate(plotdata.fig3,
                      std.error = avg.pred.diff.kinship_LB$se.fit)
colnames(plotdata.fig3)[1:3]<-c("estimate","conf.low","conf.high")

plotdata.fig3<-list(
  tidy = plotdata.fig3,
  glance = NA
)
class(plotdata.fig3)<-"modelsummary_list"

#building data for figure 3
fig3<-modelplot(list("model3"=plotdata.fig3),
                size=0.4)+
  labs(x = NULL)+
  scale_x_continuous(
    limits = c(-1, 1.25),
    breaks = seq(-1, 1.25, by = 0.25),
    labels = function(x) ifelse(x %% 0.25 == 0, x, "")) +
  geom_vline(xintercept = 0, linetype = 2)+
  coord_flip()
fig3

#saving fig3 as an TIFF file.
tiff("figure3.tiff", width = 4, height = 4, units = "in", res = 600)
print(fig3)
dev.off()

#### figure 4 ####
#extracting data for unconditional from the result(cjoint_ethnocracy)

ethnoc_unconditioned <- summary(cjoint_ethnocracy)$amce %>%
  as.data.frame() %>%
  select(Attribute, Estimate, 'Std. Err') %>%
  rename(
    term = Attribute,
    estimate = Estimate,
    std.error = 'Std. Err'
  )

ethnoc_unconditioned<-mutate(ethnoc_unconditioned, conf.low = estimate - 1.96*std.error,
                             conf.high = estimate + 1.96*std.error)
gl<-data.frame(
  Obs.="6342",
  Respondents="451"
)
ethnoc_unconditioned<-list(
  tidy = ethnoc_unconditioned,
  glance = gl
)
class(ethnoc_unconditioned)<-"modelsummary_list"

# extracting data for ethnocracy=0
ethno_out<-summary(cjoint_ethnocracy)
ethnoc_cond_0<-ethno_out$Ethnocracy1amce %>%
  as.data.frame() %>%
  select(Attribute, Estimate, 'Std. Err') %>%
  rename(
    term = Attribute,
    estimate = Estimate,
    std.error = 'Std. Err'
  )

ethnoc_cond_0<-mutate(ethnoc_cond_0, conf.low = estimate - 1.96*std.error,
                      conf.high = estimate + 1.96*std.error)

ethnoc_cond_0_list<-list(
  tidy = ethnoc_cond_0,
  glance = gl
)
class(ethnoc_cond_0_list)<-"modelsummary_list"

#extracting data for ethnocracy=1
ethnoc_cond_1<-ethno_out$Ethnocracy2amce %>%
  as.data.frame() %>%
  select(Attribute, Estimate, 'Std. Err') %>%
  rename(
    term = Attribute,
    estimate = Estimate,
    std.error = 'Std. Err'
  )

ethnoc_cond_1<-mutate(ethnoc_cond_1, conf.low = estimate - 1.96*std.error,
                      conf.high = estimate + 1.96*std.error)

ethnoc_cond_1_list<-list(
  tidy = ethnoc_cond_1,
  glance = gl
)
class(ethnoc_cond_1_list)<-"modelsummary_list"

#building data for figure 4
cm<-c("residence" = "Residence","religion" = "Religion","occupation" = "Occupation",
      "number" = "Number","nationality" = "Nationality","education" = "Education",
      "duration" = "Duration","dependents" = "Dependents","age" = "Age")
fig4<-modelplot(list("Ethnocracy=1" = ethnoc_cond_1_list,
                     "Ethnocracy=0" = ethnoc_cond_0_list,
                     "Unconditioned" = ethnoc_unconditioned),
                coef_map = cm, size=0.2)+
  geom_vline(xintercept = 0, linetype = 2) +
  labs(x = NULL)+
  scale_x_continuous(
    limits = c(-0.06, 0.125),
    breaks = seq(-0.05, 0.1, by = 0.025),
    labels = function(x) ifelse(x %% 0.05 == 0, x, "")
  ) +
  scale_color_manual(
    values = c("black", "grey40", "grey70")
  ) +
  theme(legend.position = "bottom",
        legend.title = element_blank())

fig4

#saving fig4 as a TIFF file.
tiff("figure4.tiff", width = 4.5, height = 6.5, units = "in", res = 600)
print(fig4)
dev.off()

#### figure 5 ####
#extracting unconditional from the result (cjoint_kinship) to pass to the modelplot function

kinship_unconditioned <- summary(cjoint_kinship)$amce %>%
  as.data.frame() %>%
  select(Attribute, Estimate, 'Std. Err') %>%
  rename(
    term = Attribute,
    estimate = Estimate,
    std.error = 'Std. Err'
  )

kinship_unconditioned<-mutate(kinship_unconditioned, conf.low = estimate - 1.96*std.error,
                              conf.high = estimate + 1.96*std.error)
gl<-data.frame(
  Obs.="6342",
  Respondents="451"
)
kinship_unconditioned<-list(
  tidy = kinship_unconditioned,
  glance = gl
)
class(kinship_unconditioned)<-"modelsummary_list"

#extracting data for kinship=0

kinship_out<-summary(cjoint_kinship)
kinship_cond_0<-kinship_out$Kinship1amce %>%
  as.data.frame() %>%
  select(Attribute, Estimate, 'Std. Err') %>%
  rename(
    term = Attribute,
    estimate = Estimate,
    std.error = 'Std. Err'
  )

kinship_cond_0<-mutate(kinship_cond_0, conf.low = estimate - 1.96*std.error,
                       conf.high = estimate + 1.96*std.error)

kinship_cond_0_list<-list(
  tidy = kinship_cond_0,
  glance = gl
)
class(kinship_cond_0_list)<-"modelsummary_list"

#extracting data for Kinship=1
kinship_cond_1<-kinship_out$Kinship2amce %>%
  as.data.frame() %>%
  select(Attribute, Estimate, 'Std. Err') %>%
  rename(
    term = Attribute,
    estimate = Estimate,
    std.error = 'Std. Err'
  )

kinship_cond_1<-mutate(kinship_cond_1, conf.low = estimate - 1.96*std.error,
                       conf.high = estimate + 1.96*std.error)

kinship_cond_1_list<-list(
  tidy = kinship_cond_1,
  glance = gl
)
class(kinship_cond_1_list)<-"modelsummary_list"

#drawing figure 5
cm<-c("residence" = "Residence","religion" = "Religion","occupation" = "Occupation",
      "number" = "Number","nationality" = "Nationality","education" = "Education",
      "duration" = "Duration","dependents" = "Dependents","age" = "Age")
fig5<-modelplot(list("Kinship=1" = kinship_cond_1_list,
                     "Kinship=0" = kinship_cond_0_list,
                     "Unconditioned" = kinship_unconditioned
),coef_map = cm,
size=0.2)+
  geom_vline(xintercept = 0, linetype = 2) +
  labs(x = NULL)+
  scale_x_continuous(
    limits = c(-0.0825, 0.11),
    breaks = seq(-0.05, 0.1, by = 0.025),
    labels = function(x) ifelse(x %% 0.05 == 0, x, "")
  ) +
  scale_color_manual(
    values = c("black", "grey40", "grey70")
  ) +
  theme(legend.position = "bottom",
        legend.title = element_blank())

fig5

#saving fig5 as an TIFF file.
tiff("figure5.tiff", width = 4.5, height = 6.5, units = "in", res = 600)
print(fig5)
dev.off()

#### figure 6 ####
#extracting data for unconditional from the result (cjoint_ethnocracy) to pass to the modelplot function

LB_unconditioned <- summary(cjoint_LB)$amce %>%
  as.data.frame() %>%
  select(Attribute, Estimate, 'Std. Err') %>%
  rename(
    term = Attribute,
    estimate = Estimate,
    std.error = 'Std. Err'
  )

LB_unconditioned<-mutate(LB_unconditioned, conf.low = estimate - 1.96*std.error,
                         conf.high = estimate + 1.96*std.error)
gl<-data.frame(
  Obs.="6342",
  Respondents="451"
)
LB_unconditioned<-list(
  tidy = LB_unconditioned,
  glance = gl
)
class(LB_unconditioned)<-"modelsummary_list"

# extracting data for left_behind=0

LB_out<-summary(cjoint_LB)
LB_cond_0<-LB_out$LeftBehind1amce %>%
  as.data.frame() %>%
  select(Attribute, Estimate, 'Std. Err') %>%
  rename(
    term = Attribute,
    estimate = Estimate,
    std.error = 'Std. Err'
  )

LB_cond_0<-mutate(LB_cond_0, conf.low = estimate - 1.96*std.error,
                  conf.high = estimate + 1.96*std.error)

LB_cond_0_list<-list(
  tidy = LB_cond_0,
  glance = gl
)
class(LB_cond_0_list)<-"modelsummary_list"

# extracting data for left_behind=1

LB_cond_1<-LB_out$LeftBehind2amce %>%
  as.data.frame() %>%
  select(Attribute, Estimate, 'Std. Err') %>%
  rename(
    term = Attribute,
    estimate = Estimate,
    std.error = 'Std. Err'
  )

LB_cond_1<-mutate(LB_cond_1, conf.low = estimate - 1.96*std.error,
                  conf.high = estimate + 1.96*std.error)

LB_cond_1_list<-list(
  tidy = LB_cond_1,
  glance = gl
)
class(LB_cond_1_list)<-"modelsummary_list"

# drawing figure 6
cm<-c("residence" = "Residence","religion" = "Religion","occupation" = "Occupation",
      "number" = "Number","nationality" = "Nationality","education" = "Education",
      "duration" = "Duration","dependents" = "Dependents","age" = "Age")
fig6<-modelplot(list("Left_Behind=1" = LB_cond_1_list,
                     "Left_Behind=0" = LB_cond_0_list,
                     "Unconditioned" = LB_unconditioned
),coef_map = cm,
size=0.2)+
  geom_vline(xintercept = 0, linetype = 2) +
  labs(x = NULL)+
  scale_x_continuous(
    limits = c(-0.0825, 0.12),
    breaks = seq(-0.05, 0.1, by = 0.025),
    labels = function(x) ifelse(x %% 0.05 == 0, x, "")
  ) +
  scale_color_manual(
    values = c("black", "grey40", "grey70")
  ) +
  theme(legend.position = "bottom",
        legend.title = element_blank())

fig6

#saving fig6 as a TIFF file.
tiff("figure6.tiff", width = 4.5, height = 6.5, units = "in", res = 600)
print(fig6)
dev.off()

#### end  ####
